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An  implicit  depth-averaged  2-D  finite  volume  model  has  been  developed  to  simulate  sediment  transport  and  bed 
morphological  changes  under  actions  of  currents  and  waves  near  coastal  inlets.  The  model  computes  the  depth- 
averaged  2-D  shallow  water  flow  and  non-equilibrium  transport  of  total-load  sediment,  accounting  for  the  effects  of 
wave  radiation  stresses  and  turbulent  diffusion  induced  by  currents,  waves  and  wave  breaking.  The  model  uses  a 
quadtree  rectangular  mesh  to  locally  refine  the  mesh  around  structures  of  interest  or  where  the  topography  and/or 
flow  properties  change  rapidly.  The  grid  nodes  are  numbered  by  means  of  an  unstructured  index  system  for  more 
flexibility  of  mesh  generation.  The  SIMPLEC  algorithm  is  used  to  handle  the  coupling  of  water  level  and  velocity  and 
the  Rhie  and  Chow’s  (1983)  momentum  interpolation  method  is  adopted  to  determine  the  intercell  fluxes  on  non- 
staggered  grid.  Well-developed  longshore  current  and  wave  setup  determined  with  the  reduced  1-D  momentum 
equations  are  used  as  the  cross-shore  boundary  conditions.  The  model  has  been  tested  in  several  laboratory  and  field 
cases,  showing  good  performance.  In  particular,  it  can  use  a  long  time  step  and  is  efficient  in  computation  on  a  PC 
platform.  It  has  a  potential  for  simulation  of  long-term  coastal  morphodynamic  processes. 
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INTRODUCTION 

A  coastal  inlet  connects  an  ocean  or  sea  through  a  typically  narrow  landmass  to  a  bay,  estuary, 
lagoon,  or  river.  It  is  usually  a  vital  navigation  link  and  central  for  exchange  of  water,  sediments,  and 
nutrients  between  its  two  sides.  It  often  experiences  multiple-scale  interacting  forces,  such  as  waves, 
wind,  tide,  river  flow,  and  density  current,  and  undergoes  very  complex  morphodynamic  processes. 
Prediction  of  these  processes  has  been  a  challenging  but  crucial  task  for  coastal  sediment  management 
and  navigation  channel  maintenance. 

Numerous  computational  models  have  been  developed  in  recent  decades  for  coastal  hydrodynamic 
and  morphodynamic  processes  (e.g.,  Chesher  et  al.  1993,  Roelvink  and  Banning  1994,  Ranasinghe  et  al. 
1999,  Cayocca  2001,  Fortunato  and  Olveira,  2003,  Warner  et  al.  2008).  Most  of  these  models  are  based 
on  the  assumption  that  bed  load  or  suspended  load  are  instantaneously  in  equilibrium  with  the  forcing 
processes  on  each  computational  node.  However,  because  of  the  dynamic  nature  of  currents  and  waves 
on  the  coast,  neither  bed  load  nor  suspended  load  is  usually  in  an  equilibrium  state.  The  non¬ 
equilibrium  transport  modeling  is  widely  used  in  river  sedimentation  (Han  1980,  Phillips  and 
Sutherland  1989,  and  Wu  2008)  and  has  been  recently  extended  to  coastal  sedimentation  (e.g.,  Sanchez 
and  Wu  2010).  This  approach  renounces  the  assumption  of  local  equilibrium  and  simulates  sediment 
transport  more  realistically. 

In  the  past  years,  the  Coastal  Modeling  System  (CMS)  for  nearshore  hydrodynamic  and 
morphodynamic  processes  has  been  developed  under  the  Coastal  Inlets  Research  Program  (CIRP)  of  U. 
S.  Army  Corps  of  Engineers  (Militello  et  al.  2004,  Buttolph  et  al.  2006,  Sanchez  and  Wu  2010).  The 
existing  CMS  hydrodynamic  model  (called  CMS-Flow)  solves  the  2-D  shallow  water  equations  and 
equilibrium/non-equilibrium  total-load  sediment  transport  using  an  explicit  finite  volume  method 
(FVM)  on  structured  rectangular  mesh.  The  CMS-Flow  is  two-way  coupled  with  the  spectral  wave 
transformation  model  called  CMS- Wave,  which  solves  the  steady-state  wave-action  balance  equation 
on  a  Cartesian  grid  with  a  finite  difference  scheme  and  considers  wind  wave  generation  and  growth, 
diffraction,  reflection,  dissipation  due  to  bottom  friction,  white  capping  and  breaking,  wave -wave  and 
wave-current  interactions,  wave  runup,  wave  setup,  and  wave  transmission  through  structures  (Lin  et  al. 
2008). 

To  improve  the  CMS-Flow  model’s  computational  efficiency,  an  implicit  depth-averaged  2-D 
scheme  has  been  recently  implemented  in  the  CMS  to  solve  the  shallow  water  flow  and  sediment 
transport  equations  based  on  a  quadtree  rectangular  mesh.  The  flow  equations  are  solved  within  the 
FVM  with  the  SIMPLEC  algorithm  (van  Doormaal  and  Raithby  1984)  on  non-staggered  grid  to  handle 
the  coupling  of  water  level  and  velocity,  as  presented  by  Wu  et  al.  (2010)  in  detail.  This  implicit  model 
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has  been  enhanced  by  implementing  several  new  approaches  for  model  parameters  and  capabilities, 
such  as  wave/current  interaction,  bottom  friction,  cross-shore  boundary  condition,  and  sediment 
transport.  These  enhancements  are  described  in  this  paper. 


FLOW  AND  SEDIMENT  TRANSPORT  EQUATIONS 

The  depth-averaged  shallow  water  equations  in  the  Cartesian  coordinate  system  are  written  as: 
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where  t  is  the  time;  x  and  y  are  the  horizontal  Cartesian  coordinates;  h  is  the  total  flow  depth;  u  and  v 
are  the  depth-averaged  flow  velocities  in  x-  and  y-directions,  respectively;  rj  is  the  water  surface 
elevation  above  the  reference  sea  level;  g  is  the  gravitational  acceleration;  p  is  the  density  of  flow;  vt 

is  the  eddy  viscosity  due  to  turbulence;  rSx  and  rSy  are  the  wave  radiation  stresses,  rbx  and  rby  are  the 
bed  shear  stresses,  rwx  and  rwy  are  the  wind  driving  forces,  all  in  the  x-  and  y-directions,  respectively; 
and  fc  is  the  Coriolis  force  coefficient.  Determination  of  the  wind  driving  force  and  the  wave  radiation 
stresses  refers  to  Buttolph  et  al.  (2006)  and  Lin  et  al.  (2008). 

The  moving  sediment  (total  load)  in  the  water  column  is  traditionally  divided  into  suspended  load 
and  bed  load.  The  bed  load  moves  by  rolling,  sliding  and  saltating  in  a  thin  layer  of  a  few  particle  sizes 
above  the  bed,  whereas  the  suspended  load  is  transported  by  the  turbulent  flow  in  the  water  column 
above  the  bed-load  layer.  To  reduce  the  number  of  differential  equations  to  be  solved,  the  established 
model  combines  the  bed-load  and  suspended  load  as  total  load  (bed-material  load)  and  uses  the 
following  equation  to  compute  the  transport  of  the  total  load  and  the  resulting  bed  changes: 
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where  xt  is  the  horizontal  coordinate  with  subscript  i=  1  and  2  representing  the  x-  and  y-directions, 
respectively;  ut  is  the  depth-averaged  current  velocity  in  the  z-th  direction;  Ct  is  the  depth-averaged 
total-load  concentration;  is  a  correction  factor  as  the  ratio  of  depth-averaged  sediment  and  flow 
velocities;  vs  is  the  sediment  diffusion  coefficient,  assumed  to  be  proportional  to  the  turbulent  eddy 
viscosity  as  vs  —  vt  /  crs  with  crs  being  the  Schmidt  number;  rs  is  the  ratio  of  suspended  load  to  total 
load,  which  appears  in  the  diffusion  term  to  account  for  only  the  diffusion  of  suspended  load  and  is 
determined  with  the  van  Rijn  sediment  transport  equations;  Of  is  the  sediment  fall  velocity;  C,*  is  the 

depth-averaged  total-load  concentration  at  the  equilibrium  state,  determined  with  one  of  the  sediment 
transport  formulae  available  in  the  model,  developed  by  Watanabe  (1987),  Soulsby  (1997) ,  Lund-CIRP 
(Camenen  and  Larson  2007),  or  van  Rijn  (2007a,  b);  Uc  is  the  depth-averaged  resultant  current  velocity 

v2  ;  at  is  the  total-load  adaptation  coefficient,  determined  by  at  -  U c  h  /  (Ltcof) ,  with  Lt 

being  the  total-load  adaptation  length;  f  is  the  bed  surface  elevation;  and  Ds  is  a  function  of  the  flow 
and  sediment  characteristics. 

The  depth-averaged  total-load  concentration  in  Eq.  (4)  is  defined  as  Ct  -  £  uccdz  ju ch  ,  where  c  is 

the  local  sediment  concentration  and  uc  is  the  stream- wise  local  current  speed.  With  this  definition,  the 
sediment  transport  is  simply  qt  =UchCt. 
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The  correction  factor  is  defined  as  uccdz  j {lJc  cdz^j  ,  which  accounts  for  the  time  lag 

(hysteresis)  between  flow  and  sediment  transport.  It  has  a  value  close  to  unity  for  fine  sediments,  but 
decreases  with  increasing  grain  size.  Its  evaluation  can  be  found  in  Wu  (2008). 

In  addition,  the  model  can  handle  hard  bottom  and  bottom  avalanching. 

BOTTOM  FRICTION 

Bed  shear  stress  experienced  by  the  current  in  the  case  where  both  waves  and  current  coexist  is 
determined  by 


Tbx=cfp(u2+v2+0.5U2wJl2u 

(6) 

Tby=cfp(u2+v2+0.5U2wJ,2v 

(7) 

where  Uwm  is  the  maximum  orbital  bottom  velocity  of  wave;  and  Cj-  -  gn2  / h112  ,  in  which  n  is 
Manning’s  roughness  coefficient. 

TURBULENCE  CLOSURE 

The  eddy  viscosity  may  be  affected  by  current,  waves,  and  wave  breaking  in  coastal  waters.  A 
simple  linear  combination  of  these  three  components  of  eddy  viscosity  is  used  in  this  study,  as 
described  below: 

Vt=Vt,c+Vt,w  +  Vt,b  (8) 

where  the  subscripts  c,  w,  and  b  indicate  current,  waves,  and  wave  breaking,  respectively.  The  current 
eddy  viscosity  vt  c  can  be  determined  using  several  turbulence  models,  including  the  depth-averaged 

parabolic  eddy  viscosity  model  and  the  modified  mixing  length  model.  The  modified  mixing  length 
model  is  the  combination  of  the  depth-averaged  parabolic  eddy  viscosity  model  and  the  mixing  length 
model  (Wu,  2008): 


t,c 


-y(a0w*^)  +(4 1^1) 


(9) 


where  |s|  =  |^2(fw/fx)2  +  2 (dv/dy)2  +  ( du/dy  +  dv/dx)2  J  ;  a0  is  an  empirical  coefficient,  set  as  k/6  ,  with 

k  being  the  von  Karman  constant;  lh  is  the  horizontal  mixing  length,  determined  by  lh  =  vmin(cm/z,y) , 
withy  being  the  distance  to  the  nearest  wall  and  cm  an  empirical  coefficient  between  0.3-1. 2.  Eq.  (9) 

takes  into  account  the  effects  of  bed  shear  and  horizontal  velocity  gradients  respectively  through  the 
first  and  second  terms  on  its  right-hand  side.  It  has  been  found  that  the  modified  mixing  length  model  is 
better  than  the  depth-averaged  parabolic  eddy  viscosity  model  that  accounts  for  only  the  bed  shear 
effect. 

Waves  contribute  significantly  to  lateral  mixing,  particularly  in  the  surf  zone.  The  eddy  viscosity 
due  to  waves  is  determined  using  the  Kraus  and  Larson  (1991)  formula: 

Vt,w=AumH  (10) 


where  A  is  an  empirical  coefficient  representing  the  lateral  mixing  strength,  H  is  the  wave  height,  and 
um  is  the  amplitude  of  the  horizontal  component  of  the  wave  orbital  velocity  at  the  bottom. 

Wave  breaking  also  generates  a  significant  amount  of  turbulent  eddies  in  the  surf  zone.  Its 
contribution  to  eddy  viscosity  is  determined  with 


vt,b = k 


^V1/3 


h,  or  vtb=  a. 


1/2 


(11) 


where  Db  is  the  energy  dissipation  rate  due  to  wave  breaking,  rs  is  the  wave  radiation  stress,  and  kb , 
ab  are  two  empirical  coefficients.  The  CMS  allows  application  of  both  formulae;  the  first  formula  is 
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used  by  several  researchers  (e.g.,  Camenen  and  Larson  2007),  whereas  the  second  formula  considers  the 
similarity  between  the  effects  of  bed  shear  stress  and  wave  radiation  stress  on  turbulence.  The  second 
formula  is  used  in  the  test  cases  presented  in  this  paper. 


BOUNDARY  CONDITIONS 

For  a  well-defined  problem  governed  by  Eqs.  (l)-(3),  the  flow  discharge  or  velocity  is  needed  at 
inflow  boundaries,  while  the  water  level  is  usually  given  at  outflow  boundaries  for  a  subcritical  flow  or 
at  inflow  boundaries  for  a  supercritical  flow.  In  the  context  of  coastal  waters,  the  water  surface 
elevation  is  usually  specified  at  seaward  boundary,  which  can  be  the  time  series  of  tide  levels  recorded 
at  the  boundary,  generated  using  the  harmonic  analysis  of  tidal  constituents,  or  predicted  using  a  larger- 
scale  regional  model.  Flow  discharge  is  usually  given  at  a  river  inflow  boundary,  if  there  is  one  such 
boundary. 

Along  a  cross-shore  boundary,  it  is  assumed  that  a  well-developed  longshore  current  exists.  Thus, 
the  y  (alongshore)  momentum  equation  can  be  reduced  from  Eq.  (3)  as  follows: 


vth 


dv_ 

dx 


+  ~(TSy  +Twy)~ 
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CfUcwV 
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where  Ucw=(u  +v  +0.5 Uwm)  .  Eq.  (12)  prescribes  v,  the  alongshore  component  of  current 
velocity,  which  is  used  as  the  boundary  condition  at  the  cross-shore  boundary.  The  cross-shore  (x) 
component  of  the  velocity  may  be  copied  from  internal  nodes. 

The  water  level  setup  due  to  waves  and  winds  at  the  cross-shore  boundary  can  be  determined  by 
assuming  a  zero  alongshore  gradient  of  water  level,  or  using  the  following  equation  reduced  from  the  x 
(cross-shore)  momentum  equation  (2): 


0  =  -pgh  —  +  TSx+Twx  (13) 

ax 

Near  rigid  wall  boundaries,  such  as  beaches  and  islands,  the  wall-function  approach  is  employed 
(Wu,  2008).  A  threshold  flow  depth  (a  small  value  such  as  0.02  m  in  field  cases)  is  used  to  judge  drying 
and  wetting.  If  the  flow  depth  on  a  node  is  larger  than  the  threshold  value,  this  node  is  considered  to  be 
wet;  otherwise,  this  node  is  dry.  Because  a  fully  implicit  solver  is  used  in  the  present  model,  all  the  wet 
and  dry  nodes  participate  in  the  solution.  Dry  nodes  are  assigned  a  zero  velocity.  The  wall-function 
approach  is  applied  on  the  water  edges  between  the  dry  and  wet  nodes. 

NUMERICAL  SOLUTION  ALGORITHMS 

The  present  model  uses  a  multiple-level  quadtree  rectangular  mesh,  which  can  locally  refine  the 
mesh  around  structures  or  in  high-gradient  regions  by  splitting  a  coarse  cell  into  four  child  cells  and 
using  as  many  levels  of  refinement  as  necessary,  as  shown  in  Figure  1.  This  quadtree  mesh  can  improve 
the  accuracy  of  the  model  with  a  relatively  small  increase  in  number  of  cells.  The  grid  points  are 
numbered  by  means  of  an  unstructured  index  system,  so  that  the  mesh  can  be  flexibly  generated  while  it 
has  the  merits  of  the  traditional  rectangular  mesh. 


Figure  1.  Example  of  quadtree  mesh  system 
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The  model  uses  the  non-staggered  grid  system  and  stores  primary  variables  u-,  v-velocity,  and 
water  level  on  the  same  set  of  grid  points.  The  convective  terms  in  the  momentum  equations  are 
discretized  using  one  of  the  hybrid  upwind/central  difference,  exponential  difference,  and  HLPA 
schemes.  The  HLPA  scheme  has  second  order  accuracy  in  space,  while  the  other  two  have  the  accuracy 
between  first  and  second  orders.  The  SIMPLEC  algorithm  with  under-relaxation  is  used  to  handle  the 
coupling  of  water  level  and  velocity  and  achieve  high  numerical  stability  and  efficiency.  Fluxes  at  cell 
faces  are  determined  using  Rhie  and  Chow’s  (1983)  momentum  interpolation  method,  to  avoid  spurious 
checkerboard  oscillations.  The  discretized  algebraic  equations  with  sparse,  asymmetric  matrices  of 
coefficients  are  solved  iteratively  using  the  flexible  GMRES  method  with  ILUT  preconditioning. 
Details  of  the  numerical  methods  can  be  found  in  Wu  et  al.  (2010). 

The  sediment  transport  model  solves  the  convection-diffusion  equation  with  the  same  finite  volume 
method  on  the  quadtree  rectangular  mesh  as  that  used  by  the  flow  model. 

MODEL  TESTING 

Four  cases  were  used  in  this  study  to  demonstrate  the  performance  of  the  developed  model.  They 
concern  the  longshore  current  generated  by  waves  on  a  sloped  beach,  tidal  flow  in  an  estuary,  erosion 
due  to  clear  water  inflow  in  a  wide  basin,  and  flow  and  sediment  transport  near  a  coastal  inlet.  The 
simulation  results  in  these  test  cases  are  described  below. 

Case  1:  Alongshore  Current  Generated  by  Waves  on  a  Beach 

Visser  (1982)  conducted  laboratory  experiments  of  monochromatic  waves  on  a  planar  beach  and 
collected  measurements  on  waves,  currents  and  water  levels.  Visser’ s  Case  4  was  used  to  validate  the 
developed  model.  The  bathymetry  consisted  of  a  1:10  slope  for  the  first  1-m  seaward  distance,  a  1:20 
slope  for  the  next  5-m  distance,  followed  by  a  5.9  m  distance  of  flat  bottom  to  the  wave  generator.  The 
incident  wave  had  a  height  of  0.078  m,  peak  period  of  1.02  sec  and  a  wave  angle  of  15.4°.  Case  4  was 
run  over  a  concrete  bed.  A  mesh  consisting  of  84  rows  and  147  columns  was  used  with  a  constant  grid 
resolution  of  0.15  m  in  the  longshore  direction  and  a  variable  grid  resolution  between  0.04  and  0.15  m 
in  the  cross-shore  direction.  A  zero  water  level  was  forced  at  the  offshore  boundary.  Figures  2-4 
compare  the  calculated  and  measured  longshore  currents,  water  levels  and  wave  heights,  respectively. 
The  calculated  water  level  and  wave  height  agree  well  with  measured  data,  with  root  mean  squared 
errors  (RMSE)  of  0.001  and  0.009  m,  respectively.  The  magnitude  of  the  longshore  current  is 
accurately  predicted,  but  the  location  of  the  maximum  longshore  current  is  not.  The  reason  is  that  the 
model  used  in  this  study  does  not  consider  the  effect  of  roller. 


Figure  2.  Comparison  of  measured  and  computed  longshore  currents  for  Visser  Case  4  (RMSE  =  0.051  m/s). 


Figure  3.  Comparison  of  measured  and  computed  water  levels  for  Visser  Case  4  (RMSE  =  0.001  m). 
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Figure  4.  Comparison  of  measured  and  computed  wave  heights  for  Visser  Case  4  (RMSE  =  0.009  m). 


Case  2:  Tidal  Flow  in  Gironde  Estuary 

The  Gironde  Estuary  is  located  in  southwestern  France.  It  receives  runoff  from  the  Garonne  River 
and  the  Dordogne  River  and  empties  into  the  Atlantic  Ocean,  as  shown  in  Figure  5.  The  water-surface 
width  varies  from  2  to  14  km,  and  the  flow  depth  in  the  navigation  channel  is  about  6-30  m.  The 
estuary  is  partially  mixed  and  macrotidal,  with  a  12  hour  and  25  minutes  tidal  lunar  period  and  a  tidal 
amplitude  of  1.5-5  m  at  the  mouth  (Fi  et  al.  1994).  The  simulation  domain  is  about  80  km  long,  from 
the  mouth  to  the  Garonne  River  and  the  Dordogne  River.  Because  the  domain  is  relatively  simple,  a 
uniform  mesh  with  a  grid  spacing  of  250^125  m  is  used.  The  data  measured  from  May  19  to  25,  1975 
are  used  to  validate  the  model.  The  computational  time  step  is  30  minutes.  The  Manning’s  n  is  set  as 
0.015. 

Figure  6  shows  the  calculated  flow  fields  in  flood  and  ebb  tides.  Figure  7  compares  the  measured 
and  calculated  water  levels  at  stations  lie  Verte  and  Richard,  and  Figure  8  compares  the  measured  and 
calculated  flow  velocities  at  stations  Pauillac  and  Richard.  The  measured  flow  velocities  are  1  m  under 
the  water  surface  and  1  m  above  the  river  bed.  The  amplitude  and  phase  of  both  water  level  and  flow 
velocity  are  well  predicted  by  the  numerical  model.  No  obvious  phase  difference  exists  between  the 
measured  data  and  simulated  results. 


Figure  5.  Sketch  of  Gironde  Estuary,  France 
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Figure  6.  Calculated  flow  patterns  in  flood  and  ebb  tides  at  Gironde  Estuary 


Figure  7.  Measured  and  calculated  water  levels  at  selected  stations  in  Gironde  Estuary 


Figure  8.  Measured  and  calculated  velocities  at  selected  stations  in  Gironde  Estuary 

Case  3:  Erosion  due  to  Clear  Water  Inflow  in  a  Basin 

Thuc  (1991)  carried  out  an  experimental  study  on  the  erosion  process  in  a  rectangular  basin  with  a 
loose  bed  due  to  clear  water  inflow  from  a  narrow  channel.  The  rectangular  basin  was  5  m  long  and  4  m 
wide.  It  was  connected  with  a  0.2  m  wide  and  2  m  long  channel  in  the  upstream  and  a  1.2  m  wide  and 
1.0  m  long  channel  in  the  downstream.  It  had  a  concrete  bed  covered  with  a  0.16  m-thick  layer  of  fine 
sand,  which  had  a  settling  velocity  of  0.013  m/s.  The  inflow  velocity  in  the  upstream  channel  was  0.6 
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m/s,  and  the  water  depth  at  the  outlet  was  0.15  m.  The  Manning  roughness  coefficient  in  the  basin  was 
estimated  as  0.03. 

The  computational  domain  is  shown  in  Figure  9.  The  mesh  consisted  of  62  rows  and  69  columns, 
with  finer  grid  spacing  around  the  centerline.  The  time  step  for  sediment  calculation  was  30  seconds. 
The  transport  equation  which  best  fit  the  measurements  was  the  Soulsby-van  Rijn  formula  (Soulsby 
1997).  Figure  10  shows  the  calculated  flow  pattern  in  the  domain  and  bed  change  contours  around  the 
inflow  region  at  the  center  part  of  the  basin  after  the  elapsed  time  of  4  hr.  Erosion  happens  due  to  the 
inflow  of  clear  water,  and  the  eroded  sediment  moves  downstream  and  deposits  forming  a  dune  feature. 
Figure  1 1  compares  the  measured  and  calculated  bed  changes  along  the  longitudinal  centerline  at  1  and 
4  hr.  The  calculated  erosion  and  deposition  depths  are  in  good  agreement  with  the  measured  data,  in 
particular  at  time  4  hr.  The  RMS  error  is  0.014  and  0.012  m  for  times  1  and  4  hr,  respectively. 


Inlet 


Outlet 


Figure  9.  Computational  domain  and  mesh  for  Thuc  (1991)  erosion  case 
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Figure  10.  Calculated  flow  field  and  bed  change  at  the  elapsed  time  of  4  hr  in  Thuc  (1991)  case 


Figure  11.  Calculated  and  measured  bed  changes  along  the  centerline  of  the  basin  in  Thuc  (1991)  case  at  1 
(RMSE  =  0.014  m)  and  4  hrs  (RMSE  =  0.012  m). 
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Case  4:  Flow  and  Sediment  Transport  in  Grays  Harbor 

Grays  Harbor  is  located  on  the  southwest  Washington  coast  about  45  miles  north  of  the  Columbia 
River.  The  estuary  has  a  wetted  surface  area  of  approximately  91  square  miles  at  mean  higher  high 
water  and  28  squares  miles  at  mean  lower  low  water.  The  main  input  of  fresh  water  is  from  the  Chehalis 
River.  The  3  mile  wide  entrance  has  two  convergent  rock  jetties  which  extend  from  spit  points,  as 
shown  in  Figure  12. 

The  computational  grid  consisted  of  67,000  cells  and  had  a  non-uniform  spacing  from  28  to  200  m. 
Both  the  wave  and  flow  models  used  the  same  grid.  The  spectral  waves  from  the  NOAA  buoy  46029 
were  input  at  the  model  boundaries  every  3  hours.  Wind  from  the  same  buoy  was  included  in  the  wave 
model.  The  hydrodynamic,  sediment  transport  and  morphologic  time  steps  were  set  to  15  minutes. 
Different  time  steps  between  5-30  minutes  were  tested  and  the  differences  were  negligible.  The 
Manning’s  n  coefficient  was  calibrated  as  0.018  using  field  measurements.  The  hydrodynamic  model 
was  forced  with  water  level  measurements  taken  at  Station  0.  The  bed  was  modeled  with  a  constant 
grain  size  of  0.3  mm.  The  Lund-CIRP  formula  was  used  to  calculate  the  sediment  transport  capacity. 
Finally  the  total  load  adaptation  length  was  set  to  a  constant  value  of  10  m.  The  27-day  period  from 
September  14  to  October  15  of  1999  was  calculated  with  CMS. 

Figures  13-15  compare  the  computed  and  measured  water  levels,  current  velocities  and  wave 
heights  at  selected  stations.  The  measured  data  were  obtained  from  Osborne  et  al.  (2002).  The 
agreement  between  those  calculated  results  and  measured  data  is  generally  good.  Figure  16(a)  shows 
the  computed  27-day  morphology  change  from  September  to  October  of  1999.  The  regions  of  erosional 
and  depositional  trends  are  approximately  indicated  by  polygons  and  ellipses.  Although  direct 
measurements  of  the  morphology  change  during  the  same  period  are  not  available,  the  morphology 
change  from  1993  to  2005  is  shown  in  Figure  16(b)  as  a  reference  of  the  morphological  trends.  One  can 
see  the  locations  of  erosion  and  deposition  in  the  estuary  are  generally  well  predicted.  Considering  the 
long-term  bed  changes  measured  had  been  affected  by  many  complicated  factors,  such  as  dredging  and 
jetty  breach,  it  is  difficult  to  compute  the  morphological  changes  in  this  long-term  period.  This  case 
study  application  of  the  CMS  model  will  be  further  discussion  in  a  subsequent  publication. 


Figure  12.  Topography  and  measurement  stations  at  Grays  Harbor. 
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Elapsed  time,  hr 


Figure  13.  Measured  and  calculated  tide  levels  at  Grays  Harbor,  where  MTL  =  Mean  Tide  Level  (from  top  to 
bottom  RMSE  =  0.087,  0.099,  0.219,  0.198  m). 
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Figure  14.  Measured  and  calculated  current  velocities  at  Grays  Harbor  (from  top  to  bottom  RMSE  =  0.181, 
0.176,  0.207,  0.184  0.153,  0.165  m/s). 
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Figure  15.  Measured  and  calculated  wave  heights  at  Grays  Harbor  (from  top  to  bottom  RMSE  =  0.325,  0.121, 
0.501,0.339  m). 


Figure  16.  (a)  Calculated  bed  change  from  September  to  October  of  1999;  and  (b)  Measured  bed  change  in 
1993-2005  at  Grays  Harbor.  General  trends  in  calculated  morphologic  change  are  represented. 


CONCLUSIONS 

This  paper  presents  a  depth-averaged  2-D  model  of  flow,  sediment  transport  and  bed 
morphological  changes  under  actions  of  currents  and  waves.  The  model  solves  the  depth-averaged 
shallow  water  equations  using  an  implicit  finite  volume  method  based  on  quadtree  rectangular  mesh.  It 
considers  the  effects  of  wave  radiation  stresses  and  turbulent  diffusion  induced  by  current,  waves  and 
wave  breaking.  The  quadtree  technology  locally  refines  the  mesh  around  structures  of  interest  or  where 
the  topography  and/or  flow  properties  change  rapidly,  and  adopts  an  unstructured  index  system  to 
number  the  grid  nodes  for  more  flexibility  of  mesh  generation.  The  model  uses  the  SIMPLEC  algorithm 
with  the  Rhie  and  Chow’s  (1983)  momentum  interpolation  to  handle  the  coupling  of  water  level  and 
velocity  on  non-staggered  grid.  It  specifies  the  cross-shore  boundary  conditions  by  assuming  well- 
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developed  longshore  current  and  wave  setup  which  are  determined  with  the  reduced  1-D  momentum 
equations. 

The  model  computes  the  non-equilibrium  transport  of  total-load  sediment,  accounting  for  sediment 
entrainment  induced  by  waves  and  current,  the  effect  of  bed  slope  on  bed-load  transport,  etc.  The 
sediment  transport  equation  is  solved  using  the  same  finite  volume  method  as  that  used  in  the  flow 
model. 

The  model  has  been  tested  extensively  in  many  laboratory  and  field  cases.  Presented  in  this  paper 
are  four  cases,  showing  good  performance  of  the  model  in  predicting  longshore  current  generated  by 
waves,  clear  water  erosion  in  a  wide  basin,  and  unsteady  tidal  flow  and  sediment  transport  in  estuaries. 
Because  of  the  unique  advanced  numerical  methods  used,  the  model  can  use  a  long  time  step  and  is 
efficient  in  computation  on  a  PC  platform. 

The  CMS  is  still  under  development  and  will  be  improved  further  in  the  future.  For  example,  the 
location  of  the  maximum  longshore  current  generated  by  waves  is  not  well  predicted  by  the  present 
version.  A  roller  model  is  being  implemented  to  improve  this.  The  transport  of  multiple-sized  sediment 
mixtures  is  also  considered  in  the  next  version  of  the  CMS2D  model.  Upgrades  to  the  model  are 
planned  to  incorporate  dredging  and  placement  of  dredged  sediment  during  the  simulation  period. 
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